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ABSTRACT 

We investigate the K and L band dayside emission of the hot- Jupiter HD 189733b with three 
nights of secondary eclipse data obtained with the SpeX instrument on the NASA IRTF. The 
observations for each of these three nights use equivalent instrument settings and the data from 
one of the nights has previously reported by Swain et al (2010). We describe an improved data 
analysis method that, in conjunction with the multi-night data set, allows increased spectral 
resolution (R/^175) leading to high-confidence identification of spectral features. We confirm the 
previously reported strong emission at ~3.3 jam and, by assuming a 5% vibrational temperature 
excess for methane, we show that non-LTE emission from the methane v% branch is a physically 
plausible source of this emission. We consider two possible energy sources that could power non- 
LTE emission and additional modelling is needed to obtain a detailed understanding of the physics 
of the emission mechanism. The validity of the data analysis method and the presence of strong 
3.3 /im emission is independently confirmed by simultaneous, long-slit, L band spectroscopy of 
HD 189733b and a comparison star. 

Subject headings: techniques: spectroscopic, methods: data analysis, planets and satellites: atmospheres, 
planets and satellites: individual (HD 189733b) 



1. Introduction 

The field of extrasolar planets is rapidly evolv- 
ing, both in terms of number of planets discov- 
ered and techniques employed in the character- 
isation of these distant worlds. In recent years, 
increasing attention has been directed to the de- 



tection and interpretation of spectroscopic sig- 
natures of exoplanetary atmospheres and was 
mainly pioneered using Spitzer and HST in- 



struments (eg. Agol et al.l l2010l; iBeaulieu et al. 
20081 2010: Charbonneau et al. 2002. 2005. l2008t 




Grillmair et al.l2008 


; Harrington et al.ll2006l 


Knutson et al. 


20071: ISnellen et al. 


2010b: ISwain et al.l 


2008L 
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2009al lbl: iTinetti et aDl2QQ7l l2Q10h . With the re- 
moval of spectroscopic capabilities for Spitzer at 
the end of the Spitzer cold-phase, increased ef- 
forts need to be undertaken to ensure spectro- 
scopic capabilities using ground-based observato- 
ries. As inarguably difficult as this is, various 

groups have succeeded in the detection of metal 

\ 1 I — — l 

lines and complex molecules (Bean e t al. 12010c 



Redfield et al 



Swain et ah 20 



20081 : ISnellen et all I20Q8L l2010at 
(J). In order to obtain the desired 



observations, different groups have developed dif- 
ferent techniques. These can be divided into three 
main categories: 

(1) Time-unresolved techniques: here usually 
one or more high signal-to-noise (SNR) spectra are 
taken in and out of transit and both in and out-of- 
transit spectra are differenced with the additional 
use of a telluric model. Care needs to be taken 
to not over-correct and remove the exoplanetary 
signal. 

(2) Time-resolved high-resolution: this is sensi- 
tive to very thin and strong emission lines where 
the exoplanet eclipse is followed with many con- 
secutive exposures and the emission line is identi- 
fied by the varying doppler shi ft of the planet as 
it transits ([Snellen et al.ll2010al ). 

(3) Time-resolved mid-resolution: as above, the 
exoplanetary eclipse is followed by many consec- 
utive exposures with a mid-resolution spectro- 
graph making this method sensitive to broad roto- 
vibrational transitions. The use of telluric cor- 
rections with a synthetic model is not necessary 
since we obtain a normalised lightcurve per spec- 
tral channel of which the transit depths constitute 
the spectral signatures. 



Her e, we re- analyse the original ISwain et al 
(|2010h data as well as three additional plane- 
tary eclipses observed with the IRTF/SpeX instru- 
ment. One eclipse, in particular, was obtained 
with a reference star in the slit. We used the 
time-resolved mid- resolution method pioneered by 
([Swain et al.l2010h with an improved methodology 
and data-preprocessing routine. The additional 
data in conjunction with the more advanced tech- 
niques adopted, secured results at higher spectral 
resolution and smaller error bars. Furthermore, we 
thoroughly tested our data to eliminate/quantify 
the residual telluric contamination. 



2. Observations and data reduction 

Secondary eclipse data of the hot- Jupiter 
HD 189733b were obtained on the nights of Au- 
gust 11th 2 007 ( previously been analysed by 



Swain et al.l ([20101 )). June 22nd 2009 and the 
12th of July 2009 using the SpeX instrument 



([Ravner et al.ll2003l ) on the NASA Infrared Tele- 
scope Facility (IRTF). The observations were 
timed to start approximately one to two hours 
before the secondary eclipse event, until one to 
two hours post-egress. The instrumental setup 
was not changed for these three nights. The raw 
detector frames were reduced using the standard 
SpeX dat a reduction package, SpexTool, available 
for IDL ([Cushing et al.l 12004 ). resulting in sets 
of 439, 489 and 557 individual stellar spectra for 
each secondary eclipse event respectively. The ex- 
traction was done using the aperture photometry 
setting with a two arc-second aperture. 

In addition we have analysed a fourth secondary 
eclipse of HD189733b observed on July 3rd 2010 
using the same instrument. As opposed to the 
other three nights, we observed HD189733b in the 
L-band only, with a single order, long-slit set- 
ting. The one arc minute slit allowed us to si- 
multaneously observe our target and a reference 
star with a K-band magnitude of 8.05 (2MASS 
20003818+2242065). For not saturating the target 
star, we kept the exposure time at 8 seconds and 
employed the standard ABBA nodding sequence 
throughout the night. Each AB set was differenced 
to remove the background and the final spectra 
were extracted using both a custom built routine 
and standard IRAF routines. We found both ex- 
tractions to yield the same results but the custom 
built routine performs better in terms of the final 
scatter observed. The flux received from the ref- 
erence star is on average 27 times less than that 
of the target. 

The secondary eclipses in the obtained raw 
spectra (from here onwards, 'raw' refers to the 
flat fielded, background corrected, wavelength cal- 
ibrated and extracted spectra) are dominated with 
systematic (telluric and instrumental) noise. Con- 
sequently, the spectral reduction step is followed 
by data de-noising and signal amplification steps 
as described in the following sections. 
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3. Extraction of the exoplanetary spec- 
trum 

We describe in the following subsections how 
the planetary signal was extracted from the raw 
spectra. With the nature of the observations be- 
ing a combined light (planet and stellar flux) mea- 
surement, we employ time-differential spectropho- 
tometry during the time of the secondary eclipse. 
Standard photometric calibration routines typi- 
cally achieve a ~1% level of photometric accu- 
racy, hence further de-noising is necessary to reach 
the required precision. We first removed the in- 
strument systematics in the data (data cleaning) 
and then we extracted the planetary signal in the 
cleaned data (spectral analysis). 

3.1. Data-cleaning 

To achieve the accuracy we need, a robust 
cleaning of the data is required. The cleaning pro- 
cess comprises three main steps: 1) Normalising 
the spectra, getting rid of flux offsets in the time- 
series and correcting for air mass variations. 2) 
Correcting wavelength shifts between spectra by 
re-aligning all spectra with respect to one reference 
spectrum. This step removes ~ 80% of outliers. 3) 
Filtering the timeseries of each spectral channel 
with adaptive wavelets. This step removes white 
and pink noise contributions at multiple passbands 
without damaging the un derlying data structure 
Persival fc Waldenl (l2QQ0h . 



3.1.1. Normalisation 

Firstly, we discarded the spectral information 
outside the intervals of 2.1 - 2.45/im and 2.9 - 
4.0/im to avoid the edges of the K and L pho- 
tometric bands respectively. Then, we corrected 
for airmass and instrumental effects. This was 
achieved in a two step process. We first calculated 
a theoretical airmass function, AF = exp(—b x 
airmass(t)), for each night and divided the data 
by this function. However, we found this pro- 
cedure insufficient since the baseline curvature is 
caused not only by the airmass but by other instru- 
mental effects (e.g. changing gravity vectors of the 
instrument). We hence additionally fitted a sec- 
ond order polynomial to the pre- and post-eclipse 
baseline of each timeseries and divided each single 
timeseries by the polynomial. Furthermore, we 
normalised each observed spectrum by its mean 



calculated in a given wavelength band (equation 

ED. 



p (X) - Fn ^ J A = 2J - 2.45/xm K - band 
™ ( ' ~ F n [ A = 2.9 - 4.0/oti L - band 



(1) 
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where F(X) is the flux expressed as a function 
wavelength, A, for each spectrum obtained, n. F n 

A 

and F n (X) is the normalised spectrum. In the case 
of an idealised instrument and constant airmass, 
the normalisation would be superfluous. How- 
ever, due to pixel sensitivity variations and bias 
off-sets on the CCD chip, the individual spectra 
need to be normalised to avoid frequent 'jumps' 
in the individual timeseri es. In the domain of 
the high- interfere nce limit (Pa giatakis et al.l 2007 : 
Swain et al. 2010), the astrophysical signal is pre- 
served. We investigated the effects of normalis- 
ing the spectrum over a whole wavelength band or 
smaller sub-sections of the spectrum and various 
combinations of both, but found the differences to 
be negligible. 

3.1.2. Spectra re-alignment and filtering 

After the normalisation, we constructed 2D im- 
ages with rows representing spectra of the planet- 
star system at a specific time, and columns rep- 
resenting timeseries for specific wavelengths (see 
figure UK). In figure [IK, the main sources of out- 
liers in individual timeseries, are miss-alignments 
by up to 4 pixels along the wavelength axis. We 
corrected this effect by fitting Gaussians to thin 
(FWHM ~ 5px) emission and absorption lines 
to estimate the line centres to the closest pixel. 
When the shift occurred for all the lines, the spec- 
trum was corrected with respect to a reference 
spectrum, i.e. the first spectrum in the series. 
Then cosmic rays were removed by a 2D median 
filter replacing 5<r outliers with the median of its 
surrounding 8 pixels. 

3.1.3. Wavelet de-noising 

Due to variations in detector efficiency, the cu- 
mulative flux of each spectrum depends on the ex- 



3 



act position of the spectrum on the detector (hor- 
izontal bands in figure UK), resulting in high fre- 
quency scatter in each individual timeseries. This 
effect was already attenuated by the normalisation 
step but further removal of systematic and white 
noise is required. Based on the de- n oising ap- 
proach proposed by iThatte et al.l (12010) , we have 
opted for a wavelet filtering of the individual time- 
series using the 'Wavelet Toolbox' in MATLAB. 
There are clear advantages to wavelet de-noising 
compared to simple smoothing algorithms. With 
wavelets we can specifically filter the data for high 
frequency 'spikes' and low frequency trends with- 
out affecting the astrophysical signal or losing tem- 
poral phase information. This allows for an effi- 
cient reduction of white and pink noise in the in- 
dividual timeseries. By contrast, smoothing algo- 
rithms, such as kernel regression, will impact the 
desired signal since these algor ithms smooth over 
the entire frequency spect rum iDonohol (|1995f ) and 
Persival fc WalderJ (|2000h . For a more detailed 



discussion see a ppendix and IThatte et al.l (|2010); 



Dono 



(|!98lh : 



101 (ll995h ^ |Persival fc Waldenl (|2QQUI ): Ibteml 
); ISardvl (|2000h . The use of the wavelet fil- 



et al.l (l2010h: 
2000h: LStein 



tering to each individual timeseries yielded a fac- 
tor of two improvement on the final error bars. 
The final results were generated with and with- 
out wavelet de-noising and found to be consistent 
within the respective errorbars. An example of the 
final de- noised data can be seen in figure QJ3 . 

3.2. Measuring the exoplanetary spec- 
trum 

After the data were de-noised as described in 
the previous subsection, we focused on the ex- 
traction of the planetary signal. We based our 
analys is on the approach described in IS wain et al.l 
(2010). The spectral emission features of a sec- 
ondary eclipse event are too small to be statisti- 
cally significant for an individual spectral channel. 
High signal to noise detections require a low spec- 
tral resolution, i.e. binning the data in A. This 
can be done more efficiently in the frequency do- 
main for reasons discussed below. Each timeseries 
Xi(t) (here i denotes the spectral channel) was re- 
normalised to a zero mean to minimise window- 
ing effects in the frequency domain. The discrete 
fast-Fourier transform (DFT) was computed for 
each timeseries and, depending on the final bin- 
ning, m number of Fourier-transformed timeseries 



were multiplied with each other and finally nor- 
malised by taking the geometric mean (equation 

E]). 



l/ra 



(2) 



where T is the discrete Fourier-transform and 
Xi(t) is the timeseries for the spectral channel i 
for m number of spectral channels in the Fourier 
product (m G Z + ). Since the input timeseries are 
always real and the Fourier transforms are Hermi- 
tian, we can take the n'th root of the real-part of 
the final product without loosing information. 

In the time-domain, this operation is equiva- 
lent to a consecutive convolution of X{ with X^+i, 
equation [3] 



(X % * X i+1 )[n] d = f ^X,[t]X, +1 [n - t] (3) 
t=i 

We can appreciate from equation [3] that one 
eclipsing timeseries is the weighting function of 
the other. The consecutive repetition of this pro- 
cess for all remaining (i-1) timeseries, effectively 
filters the convolved timeseries with the weighting 
function that is another timeseries. This has the 
effect of smoothing out noise components whilst 
prese rving the signal commo n to all the timeseries 



sets ((Pagiatakis et al. 2007). The final result of 



this process is the geometric mean of all timeseries. 
For an individual timeseries, the eclipse signal may 
not be statistically significant but the simultane- 
ous presence of the eclipse in all the timeseries al- 
lows us to amplify the eclipse signal to a statistical 
significance by suppressing the noise. The convo- 
lution theorem states that the Fourier transform 
of a convolution is equivalent to the dot product 
of the Fourier transforms 



T{Xi * X i+1 ) = k ® F(Xi) F{X i+1 ) (4) 

where (g) signifies multiplication in the Fourier 
space and k is a normalisation constant. This pro- 
cess is the base of the our analysis. 

3.2.1. Time- domain analysis 

Calculated the Fourier product, F[X(t)\, for i 
spectral channels, we can take the inverse of the 
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Fourier transform to obtain the filtered lightcurve 
signal. 



X{t)=F-\F[X{t)Y) 



(5) 



The lightcurves were then re-normalised by fit- 
ting a second-order polynomial to the out-of- 
transit baseline. W e modeled the fin al lightcurves 
with equation 8 of iMandel fc Agoll (l2QQ2h. using 
the sy stem parameters reported in iBakos et al. 
([20061 ), with the transit depth as the only free 
parameter left. 

As clear from the lightcurves presented in sec- 
tion [5j the systematic noise in the data is higher 
in areas of low transmissivity. Systematic noise 
increases the scatter of the obtained lightcurves 
as well as the error-bars of the final spectra and 
places a lower limit of m = 50 channels (~ 2.88nm) 
on the currently achievable spectral bin size. This 
is a noticeable improve ment compared to the orig- 
inal IS wain et al.l (|2010h analysis which reported a 
lower limit of m = 100 and 150 spectral channels 
for the K and L-bands respectively. 

3.2.2. Frequency- domain analysis 

The generated lightcurves are of high qual- 
ity and ready for accurate spectroscopic measure- 
ments. However, as previously mentioned, a cer- 
tain amount of periodic and systematic noise is 
still present in the timeseries. The noise residu- 
als are in part generated during the conversion of 
the data from the frequency domain to the time 
domain, in part are due to systematics. We can 
remove some of these residuals, by measuring the 
eclipse depth directly in the frequency domain, as- 
suming that most systematic noise is found at dif- 
ferent frequencies to the eclipse signal. 

In first order approximation, we can assume 
the eclipse signal to be a box-shaped function 
or square wave of which the F ourier transform i s 
the well known sine function (jRilev et al.ll2004f ). 
The Fourier series of such a symmetric square 
wave is given by equation Has a function of the 
lightcurve's transit depth, 5, and transit duration, 

T. 



f{t) = 4,5 ( sin(2/r) - ^sin(3/r) + Jsin(5/r) - ... 



(6) 



45 £ 



fc=l,3,5... 



sin((2fc - l)2/r) 
(2k - 1) 



The lightcurve signal is composed of a series of 
discrete frequencies, fc, since the boundary condi- 
tions of the function are finite. This series is very 
rapidly converging. Figure [2] illustrates this. Here 
we took the Fourier transform of the secondary 
eclipse model shown in the insert. The frequency 
spectrum is centred on the first Fourier coefficient. 
It is clear that most of the power is contained in 
the first Fourier coefficient and the series quickly 
converges asymptotically to zero after the third 
coefficient. Taking the product in equation [2] has 
the effect of strengthening the eclipse signal, whilst 
weakening the noise contribution: the frequencies 
contributing to the noise are in fact expected to 
be different to the ones contributing to the eclipse 
signal. In the case of stochastic (Gaussian) noise, 
wavelength dependent instrumental noise or scin- 
tillation noise, this is obvious. 

Following Fourier series properties, the mod- 
ulus of the amplitude, of the coefficients in 
equation [6] is directly proportional to the tran- 
sit depth S and the transit duration r, where 
r = ti-4/ts and ti-4 is the transit duration from 
the first to fourth contact point and t s is the sam- 
pling rate (ie. exposure time + overheads, fig. [2]). 



\A\. 
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fc=l,3,5.. 



(2fc-l) 



(7) 



The amplitude of the Fourier coefficients above 
k = 1 decreases by (2fc — 1) for a box-shape func- 
tion and is an even faster converging series for real 
lightcurve shapes which are used in the analysis 
(see appendix). Following from equation [7] we see 
that for the first Fourier coefficient, k = 1, the 
relationship between the transit depth, 5, and the 
Fourier coeffcient amplitude, \A\, is simply given 
by |Afc = i| = (t/2)S. From the analytical argu- 
ments presented above, we know that r is the 
transit duration (in units of number of observed 
spectra). We checked the consistency of the the- 
ory with the data, by calculating the value of 
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r numerically. To calculate r we produced sec- 
ondary eclipse curves with the transit duration 
and sampling; rate of the original IRTF data sets 
([Mandel fc Agolll2QQ2[ equation 8). We generated 
300 curves with transit depths (r) ranging from 
0.0001 to 0.1 and measured the corresponding am- 
plitude (|Afc = i|). Here, the derivative, d(\2A\)/dS 
gives us the value of r. We find r = 116 in-eclipse 
measurements, which agrees with the number of 
in-transit spectra obtained for the real IRTF data- 
sets. 

N spectra were obtained at a constant sampling 
interval of t s , giving us a sampling rate of R = l/t s 
in the frequency domain. For a complete repre- 
sentation of the data, the sampling rate is equal 
to the Nyquist rate, R = 2B, where B is the spec- 
tral bandwidth of the Fourier transform. The total 
number of Fourier coefficients, K, is then given by 
K = 2BN . It follows that the resolution in the 
frequency domain is determined by Af = 1/N. In 
other words, the more measurements are available 
the more Fourier coefficients can be extracted to 
describe the data and consequently the frequency 
range covered by each coefficient is smaller for a 
fixed sampling rate. 

The fact that Af is finite (Af — >> for infinitely 
sampled data-sets), means that the first Fourier 
coefficient can be contaminated by remaining noise 
signals very similar in frequency. To estimate the 
error bar on this contamination, we varied the out 
of transit (oot) length N oot by 50% and calculated 
the resulting spectrum for each Af. The error is 
then estimated as the standard deviation to the 
mean of all computed spectra. 

3.3. Application to data 

We have applied the same procedure described 
in sections 13.11 & 13.21 to the four data sets. In 
addition to the individual analysis, we also com- 
bined in the frequency domain the three data sets 
recorded with the same observational technique. 
Given that the low-frequency systematics -such 
as residual airmass function, telluric water vapour 
content, seeing, etc- are significantly different for 
each individual night, by combining multiple data 
sets, we can amplify the light curve signal and re- 
duce the systematic noise. 

To generate the final K and L-band spectra, we 
chose in equation [2] m = 100 spectral channels. 



From R S p ec tra = A C entre/AA, we get a final spec- 
tral resolution of R ~ 50. Combining all three 
data-sets together (~33 spectral channels taken 
from each observed planetary eclipse) we obtain a 
spectral resolution of ~170 and ~ 185 for the K 
and L-bands respectively. We note that the spec- 
tral resolving power for the SpeX instrument, con- 
sidering the seeing, is R ~ 800. 

4. Model 

We have simulated planetary emission spec- 
tra using line by-line radia t ive tr a nsfer models as 
described in iTinetti et aD (|2005l l2QQ6h with up- 
dated line lists at the hot temp eratures from UCL 
ExoMol and new HITEMP dBarber et all [jooi; 
lYurchenko et al.ll201ll : feothman et al.ll2009f ). Un- 
fortunately accurate line lists of methane at high 
temperatures covering the needed spectral range 
are not yet available. We combined HITRAN 2008 



(jRothman et al.l 120091), and the hig h temperature 



measurements from (|Thievin et al.l 120081 ). These 



LTE-models were fitted to the spectra presented 
in section El 

Additional to the standard LTE model, we con- 
sidered possible non-LTE models to fit the pre- 
sented data. Upper atmospheres of planetary at- 
mospheres are subject to non-LTE emissions; al- 
though negligible in most part of the near infrared 
spectrum, these emissions become dominant in the 
strongly absorbing vibration bands of molecular 
constituents, like CO2 in telluric planets and CH4 
in giant planets (and Titan). A synthetic model of 
the spectrum in the L band has been adapted from 
a model of Giant Plan ets fluorescence of CH 4 de- 
veloped for ISO/SWS dDrossart et al.lll999l ). The 
main steps involved in the radiative transfer with 
redistribution of frequency in non-LTE regime can 
be summarised as follows: 

• We first calculate the solar (stellar) flux ab- 
sorbed from all bands of CH4. Although 
classical, this part of the model can be cum- 
bersome as all the main absorption bands 
corresponding to the stella flux have to be 
(in principle) taken into account. Limita- 
tions come from the knowledge of the spec- 
troscopy of the hot bands. In this model, 
the following bands are taken into account: 
Pentad (3.3 micron) Octad (2.3 micron); 
Tetradecad (1.8 micron). An estimate of the 
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accuracy of the approximation in neglecting 
hotter bands will be give n below. Follow- 
ing: an approach given by iDovennette et al.l 
(|l998h . the spectroscopy of CH4 is simplified 
by dividing the vibrational levels in stretch- 
ing and bending modes: therefore x super- 
levels (instead of the 29 potential sub-levels 
of the molecule). It is also assumed that 
for each super-level belonging to a polyad, 
thermal equilibrium is achieved within the 
population. This assumption comes from 
the observation that intra- vibrational transi- 
tions within polyads have a higher transition 
rate than inter-vibrational transitions. 

• The population of the vibrational levels is 
then calculated within each " super- level" 
of CH4. The vibrational de-excitation is 
assumed to follow the bending; mode de- 
excitation scheme ( Applebvlll990l ). 



• From the population of the each super-level, 
the radiative rate of each level can be calcu- 
lated to determine the emission within each 
of the bands (fundamental, octad-dyad and 
tetradecad-pentad) that contributing to the 
3.3 micron domain. 

• If hot band emission can be proven to re- 
main optically thin down to deep levels of 
the atmosphere, the resonant fluorescence is 
not the same, as self-absorption is an es- 
sential ingredient of the fluorescence. Evi- 
dently, photons absorbed, on average, at a 
tau=l level have the same probability to be 
re- absorbed as re-emitted upwards. The op- 
tically thick fluorescence, including absorp- 
tion and re-emission, is therefore applied to 
the resonant band. 

5. Results 

5.1. Validation of the method used 

As described in previous sections, we analysed 
four nights of observations: three in multi-order 
mode, with only HD 189733b in the slit (re- 
ferred to as 'short-slit nights') and one night in 
L-band with single order, long-slit set up, observ- 
ing HD 189733b and a fainter reference star si- 
multaneously. While the long-slit observation cov- 
ers a narrower spectral interval compared to the 



other eclipse observations, it is a critical test of 
the methodology with its simultaneous observa- 
tions of the target and the reference star. In figure 
[TTlwe present two lightcurves: HD 189733 and the 
reference star. Both are centred at 3.31/im with 
a binning width of 50 channels (~ 2.88nm). As 
expected the HD 189733 timeseries (top) shows 
the distinctive lightcurve shape whilst the refer- 
ence star (bottom) times eries shows a null re - 
sult. We have fitted a iMandel fc Agol (|2QQ2h 
secondary eclipse lightcurve to both and found 
the HD 189733b transit depth to be 5hdi89 = 
0.0078 ± 0.0003 and S REF = 0.0 ± 0.0007 respec- 
tively. These results are in good agreement with 
the spectra presented below. 

5.2. K and L-band spectra 

The same analysis was undertaken for the three 
short-slit nights: illustrative lightcurves are pre- 
sented in figures [3] & |U In figure [3] are plotted the 
lightcurves of the 'three-nights-combined' analysis 
for the K and L-band bands centred at 2.32, 3.20, 
3.31 3.4 and 3.6 microns, with 50 channel (~ 2.88 
nm) bins. The residual systematic noise is most 
pronounced in the areas of low atmospheric trans- 
missivity, which is reflected in the error bars of 
the lightcurves and of the retrieved spectra. We 
also show the lightcurves centred on the methane 
us branch at ~3.31/im for all individual nights, 
figure [H 

Having verified the detection of HD 189733b 
eclipse in all data sets, we have generated K- and 
L-band spectra for each individual night as well 
as for all the three nights combined. The three 
individual nights are plotted in figures [5] and [3 for 
K and L bands respectively. All spectra are con- 
sistent with each other and are wi thin the error 
bars of the initial [Swain et all ( 2010f ) results. This 
said, we find the nights of the 11th of August 2007 
and July 12th 2009 of higher quality and in better 
agreement. The single night analysis supports the 
assumption that intra-night variations are negligi- 
ble which allowed us to averaged the data sets and 
hence increase the signal to noise of the final spec- 
tra. We could hence push the resolution to R ~ 
170-180 for the final combined spectra. Figures [8] 
and[9]are the three-nights-combined K and L-band 
spectra respectively. We include in these figures 
the comparison with black body emission curves 
and LTE models. It is clear from the figures that 
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the strong features observed in the L-band cannot 
be explained by standard LTE processes. 

5.3. Comparison of the observations with 
atmospheric LTE and non-LTE mod- 
els 

Even if many uncertainties subsist on the ther- 
mal vertical profile of HD189733b, the thermal 
methane emission needed to reproduce the ob- 
served spectrum would lead to brightness temper- 
atures of ~3000 K, which not only are unlikely 
given the star-planet configuration, but would also 
appear in other bands -e.g. in the vA band at 
7.8 /im- hypothesis ruled out from Spitzer obser- 
vations. While LTE models cannot explain such 
temperatures, non-LTE models with only stellar 
photons as pumping mechanism do not supply 
enough excess flux. This result is not unexpected 
since the contribution of stellar reflection from the 
planet is smaller in L band than the thermal emis- 
sion, and fluorescence is only a redistribution of 
the stellar flux (even if a small enhancement comes 
from the redistribution of frequency in the fluo- 
rescence cascade). However, a good fit can be 
obtained by assuming a vibrational temperature 
excess for methane by 5% due to an enhancement 
of the octad level population in methane which 
is higher than expected by stellar flux pumping 
(figure [9]). This increase is currently an ad- hoc 
hypothesis and simply describes the amount of vi- 
brational temperature increase required to explain 
the observed feature. 

In the case of the K-band spectrum, it is less 
obvious whether LTE or non-LTE processes are 
prevalent. We show in figure [8] a comparison with 
two LTE simulations, one including CH4 plus CO2 
in absorption as suggested by other data sets. An- 
other model was obtained with LTE emission of 
methane. However, neither of the two simulations 
perfectly capture the spectrum observed. Given 
the stronger non-LTE emission features detected 
at ~3.3 /im, one can expect to find non-LTE ef- 
fects in the K-band as well. Further observations 
are required in order to build up the required spec- 
tral resolution to decisively constrain the excita- 
tion mechanisms at work. 



6. Discussion 

In figure [5] we present the K-band spectra of 
the three separate nights. This plot shows a slight 
discrepancy between the night of the 22nd of June 
2009 compared to the other two nights analysed. 
We can observe a systematic off-set in both the 
K and L-bands (figure [7j) with this night giving 
consistently lower emission results. We associate 
this effect to the poorer observing conditions and 
a degraded quality of the data compared to the 
data obtained in the other two nights: a very high 
intrinsic scatter of the data may in fact reduce the 
eclipse depth retrieved. We estimated the average 
spectra excluding the night of June 22nd 2009 (fig- 
ure EJ) and found the results to be in good agree- 
ment with the 3 nights-combined spectrum. This 
test demonstrates the robustness of the final re- 
trieved spectrum. It should be noted that this 
issue is less severe in the L-band, since the overall 
signal strength is higher, than in the K-band. 

Whilst the K-band spectra could be explained 
with LTE models, we encounter a quite differ- 
ent picture in the L-band. The observed emission 
around ~3.3/im exhibits a very poor match with 
the predicted LTE scenario. By contrast, non-LTE 
emission of methane can capture the behaviour of 
the us branch. Similar fluorescence effects have 
been observed in our own solar system, mainly 
CO2 in telluric and CH4 in g iant solar system 
planets ([Barthlemv et al.1 12005). In section |4] we 
outline a plausible model for the creation of such 
a prominent feature. As previously mentioned, 
the increase in CH4 vibrational temperature of 
5% is presently an ad- hoc hypothesis: it simply 
describes the amount of non-LTE population re- 
quired to fit the observations, pure LTE popula- 
tions being insufficient. The source of this popu- 
lation increase can come for a variety of sources: 
XUV illumination from the star, electron precipi- 
tations, etc. which are presently not constrained at 
all. Such effects are nonetheless known in plan- 
etary physics, such as on Jupiter, where H2 vi- 
brational temperatures in the upper atmosphere 
have been demonstrated to be o ut of equilibrium 
through Ly-alpha observations ([Barthlemv et al.1 
l2005h . with a 1.4-1.5 fold increase in vibrational 
temperature. 
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6.1. Validation of observations 

The results presented here are found to be 
consistent w i th th e results initially presented by 
Swain et al.l (l2Q10h. HST/NICMOS data in the 
K-band ([ Swain et al.l l2008) and verified in the L- 
band by t he Spitzer/IRAC 3.6^m b roadband pho- 



tometry (ICharbonneau et aL 2008), see figures [6] 
& However, iMandell etHI ([201 if ), from here 
Mil, recent l y pub lished a critique of the original 
Swain et al.l (|2010h . from here S10, result report- 
ing a non-detection of any exoplanetary features 
in their analysis. Since the re sults of this publica - 
tion are in good accord with ISwain et al.l ( 2010h . 
the fundamental discrepancy between the findings 
presented here and those by Mil need to be ad- 
dressed. 

Mil argue that the L-band features reported by 
Sll were likely due to un- accounted for telluric wa- 
ter emissions rather than exoplanetary methane. 
This hypothesis poses four main questions which 
will be addressed below: (1) Do the L-band fea- 
tures look like water emissions? (2) Are the re- 
sults repeatable? (3) Do or do we not see similar 
lightcurve features in the reference star? (4) Can 
we quantify the amount of residual telluric con- 
tamination in the data? 

6.1.1. Do the L-band features look like water? 

Here the simple answer is no. As discussed in 
section [5] and shown in figures [71 & [9j the improved 
spectral resolution of these results shows that we 
are clearly dealing with methane signatures. As 
Mil pointed out, a temporary change in telluric 
opacity due to atmospheric water (or methane) 
could mimic a secondary eclipse event. However, 
for temporal atmospheric variations to mimic an 
eclipse signal in the combined result of all three 
nights, the opacity variations, as well as the air- 
mass function, would need to be identical or at 
least very similar in all data sets. The likelihood 
of such hypothesis is very small. 

In addition, we have retrieved weather record- 
ings from near-by weather stations. These include 
periodic temperature, relative humidity and pres- 
sure readings from the the cfhtQ as well as atmo- 
spheric opacity (tau) readings at 225 fim obtained 



by the CSO0 (see figure [T2]) . Spread over all three 
eclipsing events, we found no significant correla- 
tions between these parameters and the secondary 
transit shape expected. 

6.1.2. Are the results repeatable? 

A main focus throughout this publication is to 
demonstrate the repeatability of the observations. 
In section [5] we present spectra retrieved for each 
individual observing run of the three 'short-slit' 
nights and found them consistent with each other 
within the error-bars. For the methane band 
which is the most difficult to achieve measurement 
we present lightcurves for all three observing runs 
considered, figure HJ These do vary in quality from 
night to night but are found to be consistent with 
one another over a measured timescale ranging 
from August 11th 2007 to July 12th 2009. This 
test of repeatability is of paramount importance 
in asserting the validity of the analysis as a whole. 

6.1.3. Do or do we not see similar lightcurve 
features in the reference star? 

We do not see any lightcurve features in the ref- 
erence star's timeseries. As described in previous 
sections, we have obtained a fourth night in addi- 
tion to the three main nights analysed here. This 
fourth night was taken in the single-order, L-band 
only mode with a one arc- minutes long slit. This 
allowed us to simultaneously observe the target 
HD189733b and a fainter reference star, 2MASS 
20003818+2242065, over the course of a secondary 
eclipse on July 3rd 2010. We have equally applied 
the same routines outlined in section [3] to both, 
the target and the reference. In figure [T0l we plot 
the resulting lightcurves of both stars centred at 
3. 31 /mi using the standard 50 channel bin. We 
find the transit depth for HD189733b to be within 
the error bars of the other nights analysed, whilst 
the reference star timeseries is flat. Hence, the 
routines used produce a null result where a null 
detection is expected. 

Furthermore it is important to note that a 
faulty background subtraction would have much 
stronger effects on the fainter reference star than 
on the target, as any residual background is a 
proportionally larger fraction of the stellar signal. 



1 http:/ /mkwc. if a. hawaii.edu/ 



2 http:/ /www. cso.caltech.edu/ 
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We find the mean observed flux for a single expo- 
sure to be Fhdi89 ~24300e~ and Fref ~900e~ 
for the target and the reference stars respectively. 
We can now state that the observed flux is a 
sum of the stellar flux and a background con- 
tribution: F observed = F star + F back . We also 

assume that the background flux, Fb ac fc, is the 
same for both stars as they were observed simul- 
taneously on the same detector. Whatever the 
value of F ba ck may be, its relative contribution 
on the overall flux would be ~27 times higher 
for Fref than for Fhdi89- Following this ar- 
gument, if we now assume the lightcurve feature 
to be due to an inadequate background correc- 
tion (as postulated by Mil), we would expect a 
~27 times deeper lightcurve signal in the reference 
star timeseries than in HD189733b. To illustrate 
the severity of this effect, we re-plotted the time- 
series presented in figure [TOl with an additional 27 
times deeper transit than that of HD 189733b un- 
derneath. Given the flat nature of the reference 
star's timeseries though, we can confidently con- 
firm an adequate treatment of telluric and other 
backgrounds. 

6.I.4. Can we quantify the residual telluric con- 
tamination in the data? 

Using the Fourier based techniques described in 
this paper, we can quantify the remaining contri- 
bution of systematic noise and the residual telluric 
components in the spectra shown in sec. [5] As de- 
scribed in section [321 we are mapping individual 
Fourier coefficients of the lightcurve signal in the 
frequency domain. Any systematic noise or tel- 
luric contamination can therefore only contribute 
to this one frequency bin. The degree of residual 
contamination by systematics on that frequency 
bin can hence be estimated by running the routine 
described in section 13.2.21 on only out-of-transit 
and only in-transit data, i.e. removing the eclipse 
signal. Figure [13] and [14] show the planet signal 
(black) and out-of-transit and in-transit measure- 
ments of the contamination in red and green re- 
spectively. We conclude that the amplitude of the 
systematic noise and the residual telluric compo- 
nent is within the error bars of the planetary sig- 
nal. 



7. Conclusion 

In this paper we present new data on the sec- 
ondary eclipse of HD 189733b recorded with the 
SpeX instrument on the IRTF. Our data analysis 
algorithm for time-resolved, ground-based spec- 
troscopic data, is based on a thorough pre-cleaning 
of the raw data and subsequent spectral analy- 
sis using Fourier based techniques. By combin- 
ing three nights of observations, with identical set- 
tings, and a further developme nt of the data ana l- 
ysis methodology presented in IS wain et al.l ( 2Q10h . 
we could to increase the spectral resolution to R 
- 175. 

We confirm the existence of a strong feature at 
~ 3.3/im, corresponding to the methane ^3 branch, 
which cannot be explained by LTE models. Non- 
LTE processes are most likely the origin of such 
emission and we propose a plausible scheme to ex- 
plain it. 

The possibility of telluric contamination of the 
data is thoroughly tested but we demonstrate that 
the residual due to atmospheric leakage is well 
within the error-bars, both by using Fourier based 
techniques and additional observations with a ref- 
erence star in the slit. This critical test demon- 
strates the robustness of our calibration method 
and its broad applicability in the future to other 
space and ground exoplanet data. 

I.P.W. is supported by a STFC Studentship. 
We like to thank the IRTF, which is operated by 
the University of Hawaii under Cooperative Agree- 
ment no. NNX-08AE38A with the National Aero- 
nautics and Space Administration, Science Mis- 
sion Directorate, Planetary Astronomy Program. 
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Fig. 1. — Zoomed in fraction of the data prior to 
the cleaning process (A) and post cleaning (B). 
Each column is a timeseries at a specific wave- 
length and each is an individual spectrum (n) 
taken at a specific time. 
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Fig. 3. — Lightcurves of the 'three-night- 
combined' analysis for the K and L bands. 
Lightcurves are offset vertically for clarity. 
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Fig;. 2. — s ho wing; power spectrum of a 
Mandel fc Ago! (|2QQ2h model lightcurve of 
HD 189733b (inset). It can clearly be seen that 
most power of the lightcurve signal is contained 
in the first Fourier coefficient. 
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Fig. 4. — Lightcuves centred at 3.31/im with a 
bin size of 50 channels (~ 2.88 nm) for the three 
individual nights and 'three-nights-combined'. 



11 



- — Aug. 11th 2007 
June 22nd 2009 
— July 12th 2009 
— Swain et al. (201 
+ HST/NICMOS 


3) 


500K 


















V 








f 



















2.05 2.1 2.15 2.2 2.25 2.3 2.35 

Wavelength (microns) 



Fig. 5. — showing K-band planetary signal for 
the three nights separate: August 11th 2007, June 
22nd 2009 and the 12th of July 2009 in blue, red 
and green respectively. The night of June 22nd 
2009 had poor observing conditions and the data 
was significantly noisier and planetary emissions 
retrieved are systematically lower fo r this night 
in both K and L-band. Results from ISwain et al. 
( 2Q10h are shown in black. 
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Fig. 7. — showing L-band planetary signal for the 
three nights separate: August 11th 2007, June 
22nd 2009 and the 12th of July 2009 in blue, red 
and green respectively. Similar to figure [5j the 
night of June 22nd 2009 shows a systematic lower 
emission. As described previously, this may be a 
result of t he poor data quali ty of this night. Re- 
sults from IS wain et al.l (|2Q10h are shown in black. 
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Fig. 6. — showing the combined K-band plane- 
tary signal for the nights of August 11th 2007 
and July 12th 2009 only (red), excluding the 
poor data quality of the June 22nd 2009 night. 
For comparison the spectrum of all three nights 
combined (green) is overplotted. The differ- 
ence between both spectra is small and indicates 
the night of June 22nd 2009 having a small ef- 
fect o n the overa l l resu lt. Ground-based results 
fromlSwain et al.l (l2010h and HST/NICMOS data 
( Swain et al.l 2008 ). are shown in black and purple 
respectively. 
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Fig. 8. — Three night combined K band spectrum 
compared with three black body curves at 1000, 
1500, 2000 K. Furthermore two LTE models of 
CH4 in emission (turquoise) and CH4 plus CO2 
in absorption (orange). 
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Fig. 9. — Three nights combined L-band spec- 
trum. The blue discontinuous line shows a com- 
parison of the observations with the "enhanced 
fluorescent" model; non-thermal population en- 
hancement in the octad level with a 5% increase 
of vibrational temperature of CH4. Overlaid are 
black body curves at 100, 1500, 2000, 3000 K. 
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Fig. 11. — showing on the top the observed 
lightcurve of HD189733b, beneath the simultane- 
ously observed flat timeseries of the fainter refer- 
ence star. At the bottom in red is the simulated 
reference star lightcurve expected to be observed 
under the assumption that the observed signal in 
HD 189733b is due to an imperfect background 
subtraction. The flat nature of the observed refer- 
ence star lightcurve is a strong indication that the 
background subtraction was treated adequately. 
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Fig. 10. — showing the lightcurves of the long- 
slit analysis of HD189733b and the simultaneously 
observed fainter reference star beneath, centred 
at 3.31/im with the standa rd 50 channel b i nning . 
Overplotted are two fitted Mandel fc Ago] (|2002f ) 
curves for the secondary eclipse. The HD189733b 
lightcurve is in good agreement with the other re- 
sults of this paper whilst the reference star's time- 
series is noticeably flat. 
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Fig. 12. — showing from top to bottom: Temper- 
ature (deg. C, CFHT Weather station), Rel. Hu- 
midity (%,CFHT), Pressure (mb,CFHT) and op- 
tical depth, tau (225/im, CSO) for the 12nd Aug. 
2007 (blue), 22nd June (green) and 12th July 2009 
(red). The discontinuous vertical lines mark the 
secondary transit duration. 
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Fig. 13. — showing the three night combined Ri- 
band result (black), in-transit and out-of-transit 
contamination measures are plotted in blue (dash- 
dotted) and red (dashed) respectively. It can 
clearly be seen that the contamination by telluric 
components is much smaller than the planetary 
signal and the it's amplitude lies within the sig- 
nal's error bar. 




Fig. 14. — showing the three night combined L- 
band result (black), in-transit and out-of-transit 
contamination measures are plotted in blue (dash- 
dotted) and red (dashed) respectively. It can 
clearly be seen that the contamination by telluric 
components is much smaller than the planetary 
signal and the it's amplitude lies within the sig- 
nal's error bar. 
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A. Additional notes on wavelets 



As mentioned in sect ion l3.1.3| wavelet de- noising of timeseries data has several advantages: 1) wavelet de- 
composition is a non-parametric algorithm and hence does not assume prior information on the signal or noise 
properties, making it an easy to use and objective de- noising routine; 2) contrary to smoothing algorithms 
(eg. kernel regression) high and low signal frequencies are retained; 3) temporal phase information of the 
signal is preserved during the de- and re-construction of the signal. This allows for an optimal white and 
systematic noise reduction at varying frequency pass-bands. For our purposes we use a non-linear wavelet 
shrinkage by soft-thresholding of the obtained wavelet coefficients and itera t ive re const ruction of the data. 
The in tricacies of such an approach were extensively discussed by iDonoho (|l995h and Persival fc Walden 
(|2QQQh . Using the 'Wavelet Toolbox' in MATLAB, each individual timeseries underwent a 4 level wavelet 
shrinkage using "Daubechies 4" wavelets. The wavelet coefficients were estimated for each decompo sition 
step using an heuristic form of the Stein's Unbiased Risk Estimate ( SURE) for soft-thresholding (jSteinl 
Il98lh . This allows for a MINMAX coefficient estimation ([Sardvl 2000) in cases of too low signal-to-noise 
(SNR) for the SURE algorithm. After thresholding, the timeseries were reconstructed based on the obtained 
coefficients for each timeseries. 



B. Fourier analysis 

In section [3. 2. 2\ we discuss the properties of box-shaped lightcurves in the frequency domain. It is needless 
to say that this is a gross over-simplification and that the actual secondary eclipse lightcurve is more akin to 
a trapezoid (equations IB1|) rather than a square-box. In the case of a trapezoid, we can calculate the power 
to decrease by 1/k 2 for Fourier coefficients above k=l. Hence, the Fourier series for a trapezoidal shape 
converges faster (equation IB2|) . 



Wt) = Syfa)* ( *n(l/r) + _ - ... ) (Bl) 



= 8y/(2)5 jr ( sin ( fe / 4r ) + shi(3fc/4r) 

1-1 3 K V ^ 



fe=l,3,5... 

\A\trapez = — fc2 ( B2 ) 

fc=l,3,5... 

The difference between the box-car and trapezoidal shape do not affect the linear relationship between 
spectral amplitude and transit depth. We furthermore extend the argument to limb-darkened lightcurves 
that exhibit a markedly rounder morphology. These are a natural extension to the trapezoidal case and it is 
generally true that the 'rounder' the eclipse shape, the less power is contained in Fourier coefficients above 
k = 1, and hence the series are converging even faster. 
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